function [beta2, T, T_J, F, part] = SURPartiesRobust_LogGMM(Y, X, mean, RAND, n, k)


% First estimation
beta=inv(X'*X)*X'*Y;   % Get parameter estimates
U=Y-X*beta;  % Get residuals


% Second Estimation
P2=X*inv(X'*(kron(eye(n),ones(2,2)).*(U*U'))*X)*X'; %Second step projection matrix
beta2=inv(X'*P2*X)*X'*P2*Y; %Second step parameter estimates
u2=Y-X*beta2;  % Second step residuals
U2=(reshape(u2',2,n))'; % Reshape second residuals
clear U P2;

% Compute COVAR matrix
COVAR=inv((X'*X)*(inv(X'*(kron(eye(n),ones(2,2)).*(u2*u2'))*X))*(X'*X));

% Conduct Wald test on beta parameters
[F, T, T_J] = Wald_GMM(beta2, COVAR, k); 

% Estimate partial effects
bpartial=[beta2(2:k); beta2(k+2:2*k)]; % Create beta vector excluding constant
mean_l=[mean;mean;mean];    %create matrices for partial effects of binary variables
mean_h=mean_l;
mean_h(1,4)=1;
mean_h(2,5)=1;
mean_h(3,6)=1;

[part] = partial_LogGMM(beta2, bpartial, mean, mean_l, mean_h, U2, RAND, n, k);